Electron-hole puddles in the absence of charged impurities 
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It is widely believed that carrier-density inhomogeneities ( "electron- hole puddles" ) in single- layer 
graphene on a substrate like quartz are due to charged impurities located close to the graphene sheet. 
In this Rapid Communication we demonstrate by using a Kohn-Sham-Dirac density-functional 
scheme that corrugations in a real sample are sufficient to determine electron-hole puddles on length 
scales that are larger than the spatial resolution of state-of-the-art scanning tunneling microscopy. 
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Introduction. — Graphene, a single layer of carbon 
atoms arranged in a honeycomb geometry, is a two- 
dimensional (2D) system whose carriers are subject to a 
large number of scattering mechanisms affecting its trans- 
port properties in a number of intriguing ways^*^. When 
a graphene sample produced by mechanical exfoliation is 
deposited on a substrate like Si02, it displays a maxi- 
mum mobility ^ 1.0 x 10"^ — 1.5 x 10^ cm^/(Vs). The 
main scattering mechanism limiting the mobility of such 
samples is to date still unclear and the subject of a very 
intense debate^ ^. 

Martin et a/.^ were the first to demonstrate by means 
of a single-electron transistor (SET) that close to the 
charge neutrality point the carrier density distribution 
in a graphene sheet is highly inhomogeneous. Disorder- 
induced potential fluctuations break up the electron liq- 
uid into "electron- hole puddles". These flndings have 
been subsequently conflrmed by other groups^ ^ by 
means of scanning tunneling spectroscopy (STS). The 
typical STS spatial resolution is roughly two orders of 
magnitude higher than that of the SET employed in 
Ref.l4|(W ^ 150 nm). 

Due to the linear dependence of conductivity on carrier 
densit}!!, charged impurities located near the graphene 
sheet have been early on recognized as important actorJ^ 
and invoked^ to predict electron- hole puddles. Quanti- 
tative theories of carrier-density inhomogeneities taking 
into account many-body effects have also been put for- 
warcPSEH. Despite other alternatives such as frozen rip- 
pled and resonant scatterer^ff^ttSl j^ave been proposed, 
long-range Coulomb disorder is currently the most "pop- 
ular" candidate for the main scattering mechanism lim- 
iting mobility in samples on a substrate^. 

Charged-impurity scattering as the main mechanism 
of disorder has faced, however, severe experimental (and 
theoretical) difficulties. Ponomarenko et al}^ have stud- 
ied exfoliated samples deposited on various substrates 
and found a rather weak dependence of the mobility 
on the type of substrate. In particular, the authors of 
Ref. 14 have studied transport in flakes embedded in 
media with high dielectric constants, such as glycerol, 
ethanol, and water, and measured only a small increase 
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FIG. 1: (Color online) Three-dimensional plot of the cor- 
rugated graphene sample studied in this work (experimen- 
tal data are a courtesy of V. Geringei^°). The color-coding 
of the surface labels the local value of the induced car- 
rier density 5n(r) as calculated from the Kohn-Sham-Dirac 
self-consistent theory, Eqs. ([5|-(|9|. The data in this figure 
have been obtained by setting gi — 3 eV, ttee = 0.9, and 
fie ^ 2.5 X 10^^ cm~^ (see text). 



in the mobility (at temperatures above the freezing tem- 
perature of these substances). Couto et al\^ have re- 
cently reported on low-temperature transport properties 
of graphene on SrTiOs, a well-known insulator with a di- 
electric constant varying (with temperature) in the range 
3 X 10^ < esub ^ 5 X 10^. The authors of this work have 
clearly demonstrated that i) neither the carrier mobility 
nor the amplitude of the carrier-density fluctuations Sn 
are affected by the large change in the dielectric constant 
of the substrate and ii) these quantities are practically 
identical to those measured in a typical graphene sheet 
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on Si02. 

From the theoretical point of view, we will show else- 
wher^I^ that charged impurities randomly located on a 
plane (parallel to and) at an average distance d ^ 1 nm 
from the graphene sheet^^ create extremely sharp fea- 
tures in the carrier-density spatial profile, in stark con- 
trast with the smooth inhomogeneities measured using 
STS^^. Moreover, the Dirac-point mapping procedure 
exploited in Refs. [5] and |6] fails to yield trustable re- 
sults for the reconstructed carrier density at distances 
d<2 nnP. 

Motivated by this large body of literature, in this 
Rapid Communication we demonstrate that, contrary 
to the common wisdonP^, charged impurities are not 
a necessary ingredient for the existence of electron-hole 
puddles close to charge neutrality. We establish indeed 
that smooth electron-hole puddles emerge also in the 
presence of scalar and vector potentials induced by corru- 
gations only. Carrier density inhomogeneities stemming 
from ripples and corrugations have already been stud- 
ied by a few authors^^ These studies, however, have 
focussed on artificial samples whose ripples have been 
calculated by Monte Carlo or molecular dynamics sim- 
ulations. The key added value of the present work is 
twofold: i) we study a real sample using STS experi- 
mental data^^ for the height fluctuations of a graphene 
sheet on Si02; and ii) we present an approximate theory 
that allows to calculate corrugation-induced scalar and 
vector potentials from the knowledge of the STS height- 
fluctuation maps. 

From height fluctuations to scalar and vector poten- 
tials. — We analyze the 20 nm x 20 nm corrugated 
graphene sample shown in Fig.[l] The modulations in the 
height are defined by a height-corrugation profile h{r)^ 
where r = (x, y) is a 2D vector. The function h{r) is 
known experimentall}^^. Modulations in the height lead 
to stresses and to effective scalar and gauge potentials 
which couple to the orbital degrees of freedom of the 
electron gas in the sheet thereby changing the electronic 
spectrum^^. In what follows we lay down an approximate 
theory that allows us to calculate corrugation-induced 
scalar and vector potentials from the knowledge of the 
map r ^ h{r). 

We introduce the deformation ien^o^^^ Uij = Uij(r) 

as 



Uij = -{djUi + diUj + dihdjh) 



(1) 



where Ui with i = x^y are the Cartesian components of 
the 2D displacement vector u = {ux^Uy) and dx (dy) 
is a shorthand for d/dx (d/dy). In writing Eq. (fTl) we 
have neglected two non-linear terms, i.e. {diUx)(djUx) 
and {diUy){djUy)^ which are at least one order of magni- 
tude smaller that the other terms. The only non-linear 
contribution to Uij we have retained is the last term of 
Eq. ([T]), which is of the same order of magnitude of the 
first two terms in the same equation. 

The free-energy of the lattice in the presence of defor- 
mations can be written as E[u, h] = J d^r Sei[u{r)^ h{r)] 



where the elastic free-energy density per unit area £qi is 
given h:J^^ 
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Here ~ 1 eV is the bending rigidity and A = 
2.57 eV A" and u = 9.95 eV A"' are the Lame con- 
stants of graphene^ at a temperature T = 300 K (// 
has the physical significance of shear modulus). In what 
follows we neglect the first term in Eq. ([2| since this is im- 
portant only at length scales £ < {h/\u\){K/Xy^'^ ^ 1 nm 
(estimating h ^ 1 nm and |n| ^ 0.5 A). 

The equilibrium condition in the absence of external 
forces reads dkO^ik = 0, where = SE[u^ h]/5uik = 
A ^ik '^ij(^) + Uik{r) is the stress tensoi^^. Solving 
the two equilibrium equations for i = x^y allows us to 
calculate the induced in-plane displacements u{r) and 
the deformation tensor Uij(r). In Fourier transform with 
respect to r we find: 



(A + 2/i) \q\^ 2\q\^\ 



Sij 



(3) 



where J^{q) = Y.i,k ^i^kfik{q) - = 

'^QxQyfxyiq) - qlfxxiq) - qlfyyiq) and fij{q) is the 
Fourier transform of the tensor field fij{r) = 
dih{r)djh{r). 

Scalar Vi and vector V2 = Ax — iAy potentials can 
be easily calculated from the following relation^^ Vi = 



^yy 



) and V2 = g2{ux 



yy 



), where gi 



and g2 are two coupling constants. Using Eq. (|3| we find 



^1(9) 
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A + 2/X |qr|- 



A r \ X + 1^ qI- Qy . 



Ay{q) = -2g2 



\ + 2fi |g|4 
A + qxQy 



(4) 



A + 2/X \q[ 



For the coupling constant gi we use the values ^1 = 3 eV 
and gi = 20 eV^*", while ^2 = 3c/37o/4, where f3 = 
— Slog (7o)/91og(ao) ~ 2, 70 ~ 2.7 eV is the nearest- 
neighbour hopping parameter, ao ~ 1.42 A is the carbon- 
carbon distance, and c = ii/{B\/2). For the bulk modu- 
lus {B = A+/i) we use B = 12.52 eV A-^ at T = 300 kP. 
We thus find that c ^ 0.56 at this temperature. 

The real-space scalar potential Vi (r) and the two com- 
ponents of the vector potential A{r) calculated from 
Eq. Q for ^1 = 3 eV and for the sample in Fig. [l] have 
been reported in Fig. |2j Since the experimental sam- 
ple does not respect periodic boundary conditions (which 
are used in the numerical calculations below) we actually 
work with a 40 nm x 40 nm sample which has been ob- 
tained by suitably replicating the original on#^. All nu- 
merical results shown in this Rapid Communication refer 
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FIG. 2: (Color online) Left panel: color plot of the scalar potential Vi{r) (in units of meV) calculated using Eq. (|4| with 
^1=3 eV. Central panel: the x-component Ax{r) of the vector potential (in units of meV) calculated using Eq. (|4|. Right 
panel: same as in the central panel but for the ^-component Ay{r) of the vector potential. 



to the experimentally-relevant portion of the simulation 
box. 

Self- consistent Kohn-Sham-Dirac theory of the induced 
carrier density. — The external scalar Vi{r) and vec- 
tor A{r) potentials plotted in Fig. |2] and calculated 
from Eq. Q are responsible for carrier-density inhomo- 
geneities, which can be quantified by the deviation Sn{r) 
of the local density n{r) from the "background" value 
no = 2r]/Ao-\-nc. Here 2/Ao is the density of a neu- 
tral graphene sheet, = ^ 0.052 nm^ being 
the area of the unit cell in the honeycomb lattice, and 
fic is the spatially-averaged carrier density, which can 
be positive or negative and controlled by gate voltages. 
The dimensionless parameter 77 <C 1 controls the fraction 
of TT-band electrons that are described by the massless 
Dirac fermion model^. In the numerical calculations be- 
low 77 ~ 0.1. 

Since Vi{r) and A{r) change smoothly over many lat- 
tice constants, the induced density Sn{r) can be cal- 
culated^^^ by solving a single-valley (and single-spin) 
Kohn-Sham-Dirac (KSD) equation for a two-component 

spinor ^x{r) = {(p^^\r) , (p{^\r))^ : 

{a ■ [vp + A{r)] + taVKsir)} ^x{r) = ex^x{r) . (5) 

Here cr is a 2D vector constructed with the 2x2 Pauli 
matrices ai and (J2 acting in sublattice-pseudospin space, 
V = 3joao/{2h) ^ 10^ m/s is the bare Fermi velocity, 
p = —ihVr^ Icr is the 2x2 identity matrix in pseudospin 
space, and the Kohn-Sham potential, 

VK_s{r) = V^{r) + VH(r) + V^,{r) , (6) 

is the sum of the external scalar potential Vi(r), the 
Hartree potential, and the scalar exchange-correlation 
(xc) potential. 

The (classical electrostatic) Hartree potential is given 

by 

Vu{r) = [ d?r'^^^ 5n{r') , (7) 
J e\r-r'\ 

where e = (evac + Csub)/2 is an average dielectric constant, 
^vac (^sub) being the dielectric constant of the medium 



above (below) the graphene flake. For example e ~ 2.5 
for graphene placed on Si02 (the other side being exposed 
to air), while e ^ 1 for suspended graphene. 

The third term in Fks(^), "^0(^)7 is the xc potential, 
a functional of the ground-state density, which is known 
only approximately. Following Refs . [TTI and [T8l we employ 
the local-density approximation (LDA), 

T/ LDA d[n5e^^{n)] 

where fexc(^) is the excess xc energy of a homogeneous 
2 D liq uid of massless Dirac fermions with carrier density 

The ground-state density n{r) is obtained as a sum 
over the KSD spinors ^x{r): 

n(r) = iVf^[|^(-")(r)|2 + |^f)(r)|VF(eA) , (9) 

A 

where the factor A^f = 4 is due to valley and spin degen- 
eracies and uy^E) is the usual Fermi-Dirac thermal fac- 
tor. Equation (|9]) is a self-consistent closure relationship 
for the KSD equation ([5|, since the Kohn-Sham potential 
^Ks(^) is a functional of the ground-state density n{r). 

Technical details on how to solve Eqs. ([5])-([9| are dis- 
cussed at great length in Refs. llllSi 

Numerical results and discussion. — The color cod- 
ing in Fig. [l] represents the spatial map of the calcu- 
lated induced carrier density 5n{r) for a value of the 
graphene's fine-structure constant o^ee = e^ /{hve) = 0.9 
(a value commonly used value for a graphene sheet on 
a Si02 substrate). We remind the reader that has 
the physical meaning of a dimensionless coupling con- 
stant that determines the strength of electron-electron 
interactions^. A 2D color plot of dn{r) is also reported 
in Fig. |3] for the sake of clarity. In this figure we have 
presented predictions for = 3 eV (as in Fig. [T]) but 
also for gi = 20 eV. We clearly see that the carrier den- 
sity profile dn{r) breaks into electron- hole puddles with 
extensions ranging from a few nanometers to the sample 
size. Changing the value of gi from 3 eV to 20 eV leads 
merely to a change in the amplitude of carrier-density 
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fluctuations but not in the spatial pattern of electron-hole 
puddles. Since the KSD theory includes screening due to 
TT electrons, we tend to think that one should use the un- 
screened value gi ~ 20 eV to avoid a double-counting of 
screening^^. Note also the well-deflned regions of zero 
induced density, an effect that can be traced back to 
the anomalous behavior of the xc potential in systems 
of massless Dirac fermionJ^^l^. 

A more quantitative analysis than that reported in 
Fig. [l] of the degree of correlation between topographic 
out-of-plane corrugations and carrier-density inhomo- 
geneities is shown in Fig. |3] Here we plot together with 
Sn{r) contour lines of the height map h{r). From this 
flgure one infers marginal correlations between topog- 
raphy and electron- hole puddles, as already noticed in 
Refs. I18|19l for simulated ripples. More mathematically, 
the real-space scalar and vector potentials that one de- 
rives from Eq. Q are complicated functionals^^ of the 
tensor fleld fij(r)j i.e. of the height-fluctuation map 
h{r). For example, the scalar potential, is given {modulo 
a constant) by the following highly non-local expression 



Vi{r) 
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(fr'\og{\r-r'\)F{r') , (10) 



where J^{r) = X)i,j(^ij^r - ^i^j)fij{^) is the Fourier 
transform of F{q). As a consequence, carrier-density 
inhomogeneities are not correlated in a trivial fashion 
with the height map h{r). This is most transparent 
within linear-response theory in the random phase ap- 
proximatioiP^. In this limit it is possible to show that 
the induced density in response to Vi{r) for a neutral- 
on-average graphene sheet is given by 



5n{r) 



^eff 



(11) 



where the coupling constant q^^ (with physical dimen- 
sions of inverse length) is given by 
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FIG. 3: (Color online) Top panel: 
duced carrier-density profile 5n(r) 
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Fully self-consistent in- 
(in units of 10^^ cm~^) 



in the corrugated graphene sheet shown in Fig. [T] The 
data reported in this figure have been obtained by setting 
^1=3 eV, aee = 0.9, and an average carrier density 
fie ~ 2.5 X 10^^ cm~^. The thin solid lines are contour lines 
of the height map h{r). Note that there is no simple cor- 
respondence between topographic out-of-plane corrugations 
and carrier-density inhomogeneity. Bottom panel: same as in 
top panel but for gi = 20 eV. 
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(12) 



In deriving Eq. ( pTj ) we have used that the static density- 
density response function of 2D non-interacting Dirac 
fermions is Xo(^) = —Nfq/{16hv). Eqs. (11)-(12) cap- 



ture qualitatively the main features of the numerical so- 
lution of the self-consistent KSD equation even though 
they miss some important non-linear effects. Note i) the 
intriguing formal analogy between Eq. (11) and the ex- 
pression for the classical electrostatic potential in Eq. ^ 
and ii) that the coupling constant q'^^ depends on the 
screened value of ^i, gi = gi/{l + 7rA^faee/8). Moreover, 
according to Eq. (11), a reduction of the typical height 



fluctuations h by an order of magnitude, implies a sup- 
pression of the amplitude dn of density inhomogeneities 
by two orders of magnitude, in agreement with recent 
observations for graphene on h-BN^^^. 



In summary, we have shown that in a real sample 
corrugation-induced scalar and vector potentials alone 
can in principle lead to carrier-density inhomogeneities 
with length scales that are larger than the spatial resolu- 
tion of current scanning tunneling microscopeJi^. A se- 
rious comparison between experimentally-reconstructed 
carrier-density proflles and our theoretical predictions 
may lead in a near future to achieve a better under- 
standing of the main mechanism leading to electron-hole 
puddles and limiting the mobility of unsuspended sam- 
ples. While this paper focusses on graphene sheets on 
quartz, we believe that it would be very interesting to 
carry out extensive comparisons between our t heory and 
experimental data for graphene flakes on h-Bl^^PSMl^ 
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